
# get the command arguments from perl
# rsm number
# runfolder
# run number

args=(commandArgs(TRUE))

for(i in 1:length(args)){
    eval(parse(text=args[[i]]))
}	 

library(rgl)

#runfolder <- "/home/rumen/simulation1/K=1_P=1_MU=1e-06_NU=1e-05_RSM=1_B=0.02_D=0.016_NN=100_T=7300_H=300_W=300"
#run <- 2

crypts <- read.csv(paste(runfolder,"/rep",run,".crypts",sep=""),header=FALSE)
clones <- read.csv(paste(runfolder,"/rep",run,".clones",sep=""),header=FALSE)
glog <- read.csv(paste(runfolder,"/rep",run,".log",sep=""),header=FALSE)
clones <- rbind(clones,c(-1,0,0,0,0)) # dead crypts behave as having birth.rate=0
clones[clones[,1]==0,2] <- 0.02 #fix wild type birth rate

# Get all the unique sampling time points
sampling.time <- unique(glog[,3])
sampling.gens <- glog[glog[,3]%in%sampling.time & glog[,1]==0,4]
rsmdata <- read.csv("/home/rumen/cryptsim/analysis/run.simulation.coefficients.models.csv",header=FALSE)
rsm <- 1
rsmmodel <- rsmdata[rsmdata[,1]==rsm,]
repr <- as.numeric(strsplit(as.character(rsmmodel[1,3])," ")[[1]])
surv <- strsplit(as.character(rsmmodel[2,3])," ")[[1]]
muta <- strsplit(as.character(rsmmodel[3,3])," ")[[1]]
repr.reduced <- repr[repr!=0]
relative.fitness.coefficients <- NULL 
for(i in 1:length(repr.reduced)){
 combos <- combn(repr.reduced, i)
 combos <- combos + 1
 relative.fitness.coefficients <- c(relative.fitness.coefficients, apply(combos, 2, prod))
} 
reproduction.fitness.classes <- c(1.00,sort(unique(relative.fitness.coefficients)),0.0) # 1 for wild-type crypt, 0 for dead crypt

rg.ramp <- colorRampPalette(c("darkgreen", "yellow", "red","blue"))
# Make 10 shades of green to red, dont hurt reader's eyes
ncolors <- 32
col.arr <- rg.ramp(ncolors)
col.arr <- c(col.arr,"black")
max.z <- max(round(reproduction.fitness.classes*0.02,digits=4))

# PNG files to make a movie
dir.create(paste(runfolder,"/landscape.rep",run,sep=""))

for(i in 1:length(sampling.time)){
      # Get the timepoint snapshot of the grid, only clone_id
      d <- glog[glog[,3]==sampling.time[i],2]
     
      # For every clone_id in d, obtain the birth rate of that clone from .clones
      birth.rate <- clones[match(d,clones[,1]),2]
      z <- matrix(birth.rate,nrow=300,ncol=300,byrow=TRUE)
      # Add boundig box of "normal" cells
      z <- rbind(rep(0,300),z,rep(0,300))
      z <- cbind(rep(0,302),z,rep(0,302))      
      x <- 1 * (1:nrow(z))   # 2 meter spacing (S to N)
      y <- 1 * (1:ncol(z))   # 2 meter spacing (E to W)

      cola <- col.arr[match(birth.rate,round(reproduction.fitness.classes*0.02,digits=4))]
      colp <- matrix(cola,nrow=300,ncol=300,byrow=TRUE)
      colp <- rbind(colp[1,],colp,colp[nrow(colp),])
      colp <- cbind(colp[,1],colp,colp[,ncol(colp)])      
f1 <-c( -0.6793313,  0.7338294, 0.001825878  ,  0)
f2 <-c( -0.3374919, -0.3146352, 0.887188733  ,  0)
f3 <-c(  0.6516197,  0.6020787, 0.461403012  ,  0)
f4 <-c(  0.0000000,  0.0000000, 0.000000000  ,  1)
      open3d(userMatrix=rbind(f1,f2,f3,f4),zoom=1.2,windowRect=c(0,0,580,580))
      quads3d(rbind(c(0,0,0), c(300,0,0), c(300,0,150), c(0,0,150)),
      			      col="black",front="lines",back="lines")
      quads3d(rbind(c(0,0,0), c(0,300,0), c(0,300,150), c(0,0,150)),
      			      col="black",front="lines",back="lines")
      axis3d("z-+",at=c(0,75,150),
	           labels=c(0.0,round((75*max.z)/150,digits=2),
			    round(max.z,digits=2)),col="black",line=0)
      axis3d("z-+",at=c(0,75,150),labels=c(0.0,round((75*max.z)/150,digits=2),
			    round(max.z,digits=2)),col="black",line=1)
      text3d(0,0,200,text=paste("Year",sprintf("%00.1f",sampling.time[i]/365)))
      text3d(0,0,170,text=paste("Generation",sampling.gens[i]))      
      mtext3d(text="Height (H)", "x+", line = 1)
      axis3d("x+",at=c(0,300),lab=c(300,0),color="black")
      axis3d("y+",at=c(0,300),lab=c(300,0),color="black")
      mtext3d(text="Width (W)", "y+", line = 1)
      text3d(0,290,180,text="Bifurcation rate (b)")
      mtext3d(text="Bifurcation rate (b)", "z+-", line = 1)			    
      surface3d(x,y,(z/max.z)*150, color=colp, back="cull",lit=FALSE)
      rgl.snapshot(paste(runfolder,"/landscape.rep",run,"/selectiveframe",i,".png",sep=""))
      rgl.close()

      #surface3d(x,y,(z/max.z)*150, color=colp, back="cull",lit=FALSE)
      #par3d(ignoreExtent=TRUE)
      #text3d(0,0,170,text=paste("Year",sprintf("%00.1f",sampling.time[i]/365)))
      #rgl.snapshot(paste(runfolder,"/landscape.rep",run,"/frame",i,".png",sep=""))
      #rgl.close()
}
